Phase-transitions induced by easy-plane anisotropy in the classical Heisenberg 
antiferromagnet on a triangular lattice: a Monte Carlo simulation 
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We present the results of Monte Carlo simulations for the antiferromagnetic classical XXZ model 
with easy-plane exchange anisotropy on the triangular lattice, which causes frustration of the spin 
alignment. The behaviour of this system is similar to that of the antiferromagnetic XY model on 
the same lattice, showing the signature of a Berezinskii-Kosterlitz-Thouless transition, associated to 
vortex-antivortex unbinding, and of an Ising-like one due to the chirality, the latter occurring at a 
slightly higher temperature. Data for internal energy, specific heat, magnetic susceptibility, correla- 
tion length, and some properties associated with the chirality are reported in a broad temperature 
range, for lattice sizes ranging from 24x24 to 120x120; four values of the easy-plane anisotropy are 
considered. Moving from the strongest towards the weakest anisotropy (1%) the thermodynamic 
quantities tend to the isotropic model behaviour, and the two transition temperatures decrease by 
about 25% and 22%, respectively. 



I. INTRODUCTION 

The critical behaviour of classical two-dimensional 
(2D) frustrated models has raised the interest of several 
scientists in the last yftWg- Popular realization of such 
models are HeisenbergtrHa and the XYdB antiferromag- 
net on a triangular lattice, as weU as the fully frustrated 
XY model on the square latticaafl; the role of frustration 
shows up in the particular nature of the order parameter 
in the first one, and in the presence of two kinds of sym- 
metry in the other two. In the triangular-lattice mod- 
els the minimum energy configuration (say, the ground 
state) is a ferromagnetic arrangement of the spins in each 
of three sublattices, with a relative rotation of 120° be- 
tween each other. The three vertices of each lattice pla- 
quette belong to different sublattices, and it is possible 
to associate to each plaquette a vector, the chirality, that 
defines the rotation of the spin direction around the pla- 
quette. In the XY case the (staggered) chirality turns 
into a scalar order parameter whose sign distinguishes 
between two degenerate ground states. 

In this paper we consider the XXZ triangular antifer- 
romagnet (TAF), a system that shares the symmetry of 
the XY model and is a more realistic description of a spin 
system. In addition, to treat three-dimensional spins is 
a necessary step for studying the corresponding quan- 
tum system by means of thi:: mycc- quantum self-consistent 
harmonic approximationElljMy. In the XXZ TAF the 
easy-plane anisotropy results in a double degeneracy of 
the ground state. This model is thus expected to belong 
to the universality class of the frustrated XY model. A 



very fascinating aspect of the thermodynamics of such 
a system is that it has two order parameters, the (in- 
plane) magnetization and the chirality, with two distinct 
symmetry groups, continuous 5*0(2) rotations and dis- 
crete Z2 lattice reflections, respectively. This has supKis.- 
ing consequences in view of Mermin- Wagner's theoremli3, 
that can be applied only in the first case to infer that 
the magnetization must vanish at any nonzero temper- 
ature, according-Lo the expected Berezinskii-Kosterlitz- 
Thouless (BKT)Ej critical behaviour associated with the 
rotation symmetry in the xy plane, while long-range or- 
der and an Ising-like phase transition are allowed for 
the chirality. This rich structure was observed in Monte 
Carlo (MC) simulations for the XY model. However, it 
was not clear whether the two transitions are distinct, 
with an intermediate phase, or they are manifestations 
of a new universality class, in which the two transition 
temperatures could be molten in a single, multicritical 
point: from early numerical simulationsQiJ they turned 
out to be weakly different, but in view of their uncer- 
tainty no firm conclusions could be inferred. Very re- 
cently, high-precision MC studies of the XY TAF jpodel 
and of its Villain version were reported by OlssonQ, who 
established that the BKT transition occurs at a temper- 
ature ~ 1-1.4% lower than the Ising-like one. 

Qualitatively, the results we find for the XXZ TAF 
are rather similar to those already known for the XY 
TAF, for any anisotropy strength considered, even close 
to the isotropic limit. The observed changes refiect the 
enhanced spin fluctuations out from the easy plane and 
the crossover towards the isotropic limit. The features of 
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both the Ising-hke and the BKT transition shift to lower 
temperatures. Although the numerical uncertainty is not 
much smaller than their difference, from our simulations 
the BKT critical temperature turns out slightly smaller 
than the other one, systematically for all anisotropy val- 
ues. 

In Sec. H we introduce the XXZ model and its ground 
state, discussing the connections with the correspond- 
ing XY model and its critical behaviour. 



Ill the 



In Sec. 

MC simulation algorithm is described and the definitions 
of the calculated thermodynamic quantities are given. 
Eventually, in Sec. IV the MC results are reported, gen- 
erally for four different values of the anisotropy constant, 
and analyzed for their critical behaviour and the finite- 
size effects. Conclusions are briefly drawn in Sec. 0. 



II. THE XXZ MODEL 

A. Definition of the model 

We consider the classical XXZ Hamiltonian 
J 
2 



n 



(1) 



where J is the positive (antiferomagnetic) exchange con- 
stant, and (sf,sf,sf) are the Cartesian components of 
unitary vectors, the classical spins, sitting on the sites 
{1} of a two-dimensional triangular lattice. The interac- 
tion is restricted to nearest-neighbours and d runs over 
their relative displacements (|d| = 1). The use of unitary 
spin vectors produces no loss of generality in the study 
of a classical system since the case of generic spin S can 
be taken into account by rescaling the exchange constant 
J JS^. The exchange constant sets the natural energy 
scale of the system: hence, in the following, energies and 
temperatures will be always given in units of J. 

The planar character of the system is due to the pres- 
ence of the constant A € [0, 1), which weakens the inter- 
action of the z spin components, energetically favouring 
configurations with the spins lying in the xy plane feff y- 
plane). For A = 1 the isotropic Heisenberg modeld'EI is 
recovered; for A = (when the model is also known as 
XXO) the spin components on the z axis do not even ap- 
pear in the Hamiltonian, making it formally identical to 
that of the XY.jnodel, which was the subject of some 
previous worksQB. However, in the XY and XXO models 
the phase-spaces are different: in the former the spins are 
two-dimensional vectors, while in the latter they can fluc- 
tuate out of the xy plane. Of course the thermodynamic 
properties are expected to be quantitatively different in 
the two models. 



B. Ground state properties 

The ground state configuration of can be foundEl 
by minimizing the energy of any single elementary tri- 
angular cell of the lattice. In this way one gets that. 



for every value of A G [0,1], the ground state consists 
of coplanar spins forming ±2tt/3 angles between nearest- 
neighbours (see Fig. 0). In contrast to the isotropic case, 
where the plane in which the 2tt/3 structure lies can take 
any direction in spin space, in the XXZ model (as well as 
in the XY model) such structure must take place in the 
easy-plane. In any case this leads to a VS x -v/S periodic 
ground state. 

In the XY and XXZ unfrustrated systems, like the fer- 
romagnet on a general lattice or the antiferromagnet on 
a bipartite lattice, the degeneracy of the ground state, 
connected to the possible equivalent choices of the direc- 
tion of alignment of spins in the space, corresponds to 
the S0{2) symmetry group. In the planar TAF the frus- 
tration effect causes an additional two-fold degeneracy of 
the ground state, which is due to chirality (or helicity), 
defined as the sign of rotation of the spins along the sides 
of each elementary triangle. Since global spin rotations 
in the easy-plane conserve the chirality, one configuration 
cannot be obtained from the other one by pure rotations 
but it is necessary to include some other symmetry oper- 
ation such as lattice reflection; i.e. the whole degeneracy 
corresponds to the group 50(2) x Z2. 




FIG. 1. Two degenerate ground states. The plus and minus 
denote the sign of the chirality of each elementary cell and 
the two configurations are characterized by opposite staggered 
chirality. The letters A, B and C label the three sublattices 
on which the spins are ferromagnetically aligned. 



C. Phase transitions 

As it has been pointed out by Lee et a/.i, the exis- 
tence of the extra degeneracy of the ground state allows 
the 2D frustrated XY model to support, in addition to 
spin waves and vortices, a third type of elementary ex- 
citation associated with domain-wall formation between 
regions with opposite staggered chirality (solitons). The 
first two are responsible for the loss of orientational order 
with increasing temperature, while the latter causes the 
destruction of chirality order. Then, the thermodynam- 
ics is characterized by the interplay of these three types 
of elementary excitations. 

The ground state of the three sublattices consists of 
ferromagnetically aligned spins, the interaction between 
those of the same sublattice being mediated by the 
spins of the other two. Thus the physics of the three 
sublattices is expected to be similar to the physics of 
the whole lattice in the XY unfrustrated models and a 
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Berezinskii-Kosterlitz-Tbauless (BKT) topological tran- 
sition becomes plausiblal3. Besides, since lattice reflec- 
tion is a discrete symmetry, the system can have long 
range chirality order at finite ternperature without violat- 
ing the Mermin- Wagner theoremllj and an order-disorder 
Ising-like transition is expected. Indeed this is the situ- 
ation observed for the first time, in the MC study per- 
formed by Miyashita and Shibaj for the 36-state clock 
model and, later, by Lee et ala, on the basis of group- 
theoretical symmetry arguments, combined with more 
extensive MC simulations including somewhat larger lat- 
tices (up to 72 X 72) and without the limitations intro- 
duced by the discrete sampling of the spin orientation. 
Perhaps the most striking characteristic of the phase 
transition is the logarithmic divergence of the specific 
heat, quite similar to that of the Ising model, and for 
this reason it was associated with the order-disorder tran- 
sition of the chirality. The sharp drop with increasing 
temperature of the staggered chirality, the counterpart 
of Ising magnetization, and the divergence of the corre- 
sponding chiral susceptibility are the other main char- 
acteristics of the transition connected with the loss of 
chirality order. Instead, the main sign of the vortex- 
antivortex unbinding driven BKT transition is the diyei;- 
gence of the spin correlation length and susceptibilityclfl. 
Whether the two transitions are very close but distinct, 
each conserving the proper critical exponents, or they 
melt in a single phase-transition which has both Ising 
and BKT character in a possible new universality class, 
is an open problem also shared by the other realizations 
of the so-called fully frustrated XY model, whichpare gen- 
erally believed to have . siv wlar critical behaviouiCI. 

It is well knownyyiiiU that the unfrustrated XXZ 
model shows the BKT transition, as the XY model, at 
finite temperature for every value of the anisotropy con- 
stant A < 1 (no matter how close to 1); in particular, 
for A ^ 1, the critical temperature is predicted to ap- 
proach as 1/1 ln(l — A)|. The phase transition in the 
XXZ model is again due to vortex unbinding as in the 
pure XY model, even if the vortices may have different 
character since the spins are no longer obliged to lie in 
the xy plane, with the spins in the voijtex core pointing 
preferably in the out of plane directioi£3. Likewise, it is 
at least plausible that the order-disorder transition of the 
chirality is also present in the frustrated XXZ model for 
A < 1; i.e., as far as the system shows a planar character 
and its ground state (the 27r/3 structure of the ground 
state being forced to lie in the xy plane) has the two-fold 
additional degeneracy. 

III. MONTE CARLO SIMULATION 

A. Simulation algorithm and procedure 

We performed standard MC simulations on triangu- 
lar lattices of size L x L (along the primitive vectors di- 
rections), containing N — L'^ sites and 2N elementary 
cells, with periodic boundary conditions. L was between 



24 and 120: multiples of 3 have been chosen for L to 
preserve the ground-state translation symmetry. In or- 
der to reduce the MC correlation time, a combinati on , o f 
Metropolises and over-relaxed algorithm was usedEllE3, 
as in ptpvious studies of the XXZ model on the square 
latticelia. For any given lattice size and for any temper- 
ature, the simulation procedure was the following. An 
initial configuration was generated at random and then a 
given number, typically 10 000, of Metropolis moves (one 
move consists of N single-particle moves) were made in 
order to bring the system to a thermalized configuration; 
after that, the accumulation of the averages started. A 
sample configuration, used to update the averages of the 
thermodynamic quantities, was taken after iVo — 4 over- 
relaxed moves and A'm — 2 Metropolis moves, the values 
of A^o Eind A^M being a compromise between computa- 
tional effort and minimization of MC correlation time. 
Typical runs sampled 30000 configurations. Since to 
our knowledge, no data is available for the XXZ TAF 
model, the simulation code we developed was checked in 
the isotropic case, for which we got results in leopiplete 
agreement with those reported in the literatureS'El. 

The strategy adopted to choose the lattice sizes and the 
temperatures during the simulations was the following. 
We performed a sequence of simulations for the smallest 
lattice size {L — 24) with a big step in temperature, in 
order to get a clue about the qualitative behaviour of the 
various quantities and to locate, at least approximatively, 
the critical region. Then it has been necessary to perform 
various simulations in the neighbourhood of the critical 
temperature and for larger lattice sizes, in order to locate 
the transition temperatures and to evaluate the finite- 
size effects which strongly affect the behaviour of some 
thermodynamic quantities. 

In the following, the symbol (Q) will denote the MC 
average of a general quantity Q, defined as 

M 

iQ) = MT.Q^' (2) 

M being the total number of configurations j sampled 
during the simulation. The uncertainties reported in the 
following are statistical errors, estimated in the standard 
way from the quadratic fluctuations of the correspond- 
ing observable. The effects of correlations were included 
multiplying the pure statistical errors by •\/2r, r being 
the correlation time deduced from the analysis lOf the se- 
quence of the sampled data for that observablec3. 

B. Thermodynamic observables 

We evaluated several thermodynamic quantities in or- 
der to observe the effects of the destruction of chirality 
order, such as internal energy, speciflc heat, chirality and 
its susceptibility. On the other hand, spin correlation 
functions, correlation length and susceptibility, were cal- 
culated to investigate the presence of a BKT transition. 
The internal energy per spin is defined as 
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N 



(3) 



where E = H/J. The specific heat can be computed, 
using the relation 



1 (E') {Ef 
N t2 



(4) 



with the reduced temperature t — T/ J. 

The definition of the chirahty, that is the sign of ro- 
tation of spins on the elementary triangler^H, the ground 
state configuration, is usually generalizccP-a to nonzero 
temperatures as follows: 
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(si X S2 + S2 X S3 + S3 X Si) 



(5) 



where 1,2 and 3 are the corner sites of the elementary cell 
centered at the dual-lattice vector r, always (for both up- 
ward and downward triangles) ordered in the same way 
(for example, counterclockwise). Unlike the XY model 
where, being the spins confined in the xy plane, it has 
only the component along the z axis, in the XXZ model 
the chirality is a true vector. It is normalized to 1 for 
a complete 27r/3 structure, when it has only the z com- 
ponent that can take the values ±1, as the Ising spin. 
At any finite temperature the length of compared to 
that of the other two components gives a measure of the 
rigidity of the 27r/3 structure, and can be taken as the 
order parameter. During the simulations we calculated 
the staggered chirality, defined as 
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(6) 



where the factor (— )"" assumes the values ±1 for down- 
ward/upward triangles, respectively. 

Also interesting is the susceptibility associated to stag- 
gered chirality along the z axis, calculated as 
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i-y< 



On the other hand, the sublattice in-plane spin sus- 
ceptibility was computed as the average of the sublattice 
squared magnetization, 



A=A,B,C 



iSA 



(8) 



since the average magnetization in the thermodynamic 
limit is 0. Such a definition retains its value in investi- 
gating the divergence of x ^-Iso for finite lattice simula- 
tions, because, even if different from 0, the missing term 
in Eq. (^) is well behaved and negligible in the critical 
region. Alternatively, the same information can be also 
obtained from the total k-dependcnt susceptibility. 



X(k) 



(9) 



taken at the ordering wavevector K, i.e., one of the six 
vectors pointing towards the corners of the first Brillouin 
zone of the whole lattice; for example K = (47r/3,0). A 
straightforward calculation shows indeed that the sublat- 
tice susceptibility, defined in Eq. (||), and the total one 
satisfy the following relation 



x(K)-^X-^x(k=0) 



(10) 



where x(k=0) is at t = 0, and is small with respect to 
the first, also near to the critical point. 

The spin correlation length is defined assuming the 
asymptotic exponential decay form of the in-plane spin 
correlation functions. 



C(n) = (s?s^ 



i+n 



(11) 



with i and i -I- n belonging to the same sublattice, and 
large values of n. Even if a direct fit of the two-point 
correlation function ( [ll| ) can be used, we adopted a faster 
and more reliable method to evaluate ^. This can be 
achieved translating Eq. ( |ll| ) in the reciprocal space. In 
fact, the asymptotic exponential decay in real space is 
associated with the Ornstein-Zernicke form of the Fourier 
transform of the (full lattice) spin correlation function, 
i.e., the k-dependent susceptibihty x(k), behaving as 



x(K+k) 



1 



(12) 



for small values of the wave vector k. Since the first 
Brillouin zone of the finite lattice is discrete, it is not 
possible to take arbitrarily small values of k: we used a 
fit with the first four shells around k = 0. An alternative 
way we used was to extract the value of using just the 
smallest available k: 



X(K) 
x(K+ki) 



-,1/2 



(13) 



(7) with, e.g., ki = (0,47r/i\/3). 



IV. RESULTS AND COMMENTS 
A. Thermodynamic behaviour 

We performed simulations on the XXZ model for values 
of A = 0, 0.5, 0.9, 0.99, ranging from the strongest easy- 
plane anisotropic case to the quasi-Heisenberg case. The 
isotropic model was also considered not only to compare 
the results of ous gknulation code with the data reported 
in the literatureEl'aa but also to check the consistency of 
the quasi-Heisenberg limit. 

The internal energy per spin e for the XXZ model with 
A = 0, is reported in Fig. g as a function of the reduced 
temperature t. For comparison, data for the XY model 
taken from Ref. || and for the isotropic model are also 
shown. It is evident from the figure that the qualitative 
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behaviour of the internal energy is quite different in the 
isotropic and planar cases. In fact in the latter cases the 
internal energy presents a narrow region in which size- 
dependence is apparent, and the slope becomes steeper 
the larger the lattice size. Instead, the behaviour of the 
same quantity in the isotropic case is smoother and only 
with a weak size dependence (not shown in the figure). In 
the low temperature region the internal energy, accord- 
ing to the spin wave approximation and the equipartition 
theorem, is linear in t. While in the XY model the inter- 
nal energy starts from the ground state value Cq = —3/2 
with slope 1/2, in the other two cases the slope is 1, be- 
cause of the different number of degrees of freedom: one 
per spin in the former, two in the latter cases. Increasing 
the temperature, the excitation of the out-of-plane com- 
ponent of the spins becomes more and more important 
and causes the different behaviour between the isotropic 
and XXO case. Similar behaviour is observed for all the 
values of A 7^ considered. 
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FIG. 2. Internal energy for the XXO model for lattice sizes 
L = 24 (circles), L = 36 (diamonds), L = 48 (down triangles), 
L — 60 (up triangles) and L = 120 (squares). Data for the 
XY model with L — 30, taken from Ref. ^ (full squares), and 
for the isotropic case with L = 60 (crosses) are also shown 
for comparison. In the lower figure the data are reported in 
a magnified scale in order to emphasize the size-dependence. 
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FIG. 3. Specific heat for the XXZ model for the reported 
values of A and different values of L. In the top graph the 
specific heat observed in the XY model, for L — 72, taken 
from Ref. ^ is also shown (crosses). Open symbols as in 
Fig. I 

The specific heat data are reported in Fig. || for all the 
values of A (including the isotropic case). For A 7^ 1 the 
specific heat shows the signature of a divergence which is 
an important feature of the frustrated planar antiferro- 
magnet, also present in the XY case (the corresponding 
peak, taken from Ref. is shown in the figure on the 
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top). As A ^ 1 the peak and the size dependence of 
its height become less and less pronounced, until, in the 
isotropic limit, no divergence at all is observed. The size 
dependence of the peak height is shown in Fig. ^, for 
A = 0, 0.9 and 0.99. It suggests a logarithmic divccg£;nce 
with L, just as in the two-dimensional Ising modelEJ. 
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FIG. 4. Maximum of the specific heat as a function of L 
for A = (circles), 0.9 (diamonds) and 0.99 (triangles). The 
straight lines are guides for the eye. 
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FIG. 5. Staggered chirality (a) and chiral susceptibility (b) 
for A = and various values of L, as functions of temperature. 
Symbols as in Fig. ^. 



The critical behaviour associated with the order-dis- 
order transition can be also observed in the staggered 
chirality and in its susceptibility, defined in Eqs. (^, (|^), 
and (^, and shown in Fig. |^ for A = 0. At quite low tem- 
peratures the system displays chirality order, witnessed 
by the high value of k. As the temperature raises, the 
number of cells with small chirality increases, and do- 
mains with opposite staggered chirality develop in the 
lattice. This leads to a sharp drop of the chirality and to 
the divergence of the chiral susceptibility at tc — 0.41 . 

The transition is also shown by the behaviour of the 
correlation function of the staggered chirality, namely 



(14) 



R being one of the Bravais vectors of the dual lattice. 
This is shown in Fig. |^ for A = 0, L = 60, and several 
temperatures. A low temperature phase in which the 
system is completely correlated is evident for t < 0.42. 
Increasing the temperature, the correlation functions fall 
off exponentially. Of course it is hard to extract from the 
correlation functions an accurate estimate of the criti- 
cal temperature since, approaching tc from above, when 
the correlation length becomes larger than the sampled 
lattice size, the system behaves as if it were correlated, 
eveir if, in the thermodynamic limit, it may be not. How- 
ever, we can approximatively locate the critical region, 
for A = 0, between t — 0.4 and t = 0.42, consistently 
with the behaviour of the specific heat, the chirality and 
its susceptibility. 
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FIG. 6. Chirality correlation function for A = 0, L = 60, 
and various temperatures: t = 0.3 (full circles), t — 0.4 (full 
triangles), t = 0.42 (full diamonds), t = 0.44 (circles), t = 0.5 
(triangles), t = 0.6 (diamonds). 

Let us turn now to the rotational degrees of freedom. 
We recall that the peculiarities j-ef the BKT transition 
can be summarized as foUowstd'ta. First of all, accord- 
ing to Mermin- Wagner's theorem, the system does not 
show any finite magnetization in absence of an applied 
magnetic field, at any t ^ below and above the transi- 
tion temperature ^bkt- Although the system cannot dis- 
play long-range order, below the critical temperature it 
is characterized by quasi long-range order, whose macro- 
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scopic consequence is the power-law decay of the corre- 
lation functions of the in-plane spin components, 



0.20 



1 

71'' 



(15) 



where the critical exponent ?7 is a function of temper- 
ature and assumes the universal value of 77 = 1/4 at 
^BKT- This quasi long-range order is not destroyed by the 
excitation of vortex-antivortex pairs until when, raising 
the temperature, the pairs unbind and the system under- 
goes a transition to a disordered phase with exponentially 
decaying correlation functions. The in-plane correlation 
length and susceptibility both diverge exponentially for 
t 



^ oc e 
X oc e 



{'{(t-tBKT)"^^^ 



(16) 
(17) 



(where x(K) can also be taken for x) ^-nd are infinite 
for t < tBKT- Thesa-|giroperties have been also ob- 
served in the XY TAFq'EI, which shares a similar low- 
temperature phase with the ferromagnetic counterpart, 
if the sublattice magnetization of the former replaces the 
uniform magnetization of the latter; for example, the low- 
temperature phase is well described by Eq. (15), where i 
and i -I- n belong to the same sublattice. 

Another important property of the BKT transition in 
unfrustrated planar systems is the behaviour of the spe- 
cific heat, which displays a maximum slightly above rtte 
transition temperature (usually at i ~ 1.1 1.2 ^ekt"). 
However this nuwfimum cannot be observed in frustrated 
planar systemgj'El, where it is hidden by the divergence 
connected with the chirality transition. 
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FIG. 7. Spin correlation function for A — 0.9, L = 120 and 
different temperatures: t — 0.32 (full triangles), t = 0.35 (full 
diamonds), t = 0.36 (circles), and t = 0.40 (squares). Lines 
are best fits against Eq. (|l|) (dashed lines) and Eq. (y) (full 
lines) . All fitting functions were properly symmetrized to take 
into account the periodic boundary conditions applied to the 
simulated finite-size lattice. 
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^'^(t) < 1/0.0256 ~ 39. Indeed, this is not the case for 
t — 0.35, when ^ ~ 90 and the linear fit fails. The dashed line 
is a guide for the eye. 
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Fig. displays the correlation functions of the in-plane 
spin components C(n) for A = 0.9 and various tempera- 
tures, for L = 120. Fittings against Eq. (|l|) and Eq. (0) 
are also shown in the figure. The data for t = 0.32 and 
t = 0.35 can be fitted only by the power law, while for 
t — 0.36 and t — 0.4 the exponential decay fits best: 
given the rather large lattice size, it can be reasonably 
argued that the critical temperature tBKT(A=0.9) is lo- 
cated in between between t = 0.35 and t — 0.36. As 
already noticed, it is difficult to extract accurate infor- 
mations about tBKT from the correlation functions since 
it is not possible to discriminate between the high- and 
low-temperature predicted behaviour unless data for lat- 
tice sizes L > ^ are available, requirement which cannot 
be achieved close to the critical temperature. 

As already said the methods we used to extract the 
value of the correlation length ^ were the fit, according 
to Eq. , of the total in-plane k-dependent susceptibil- 
ity x(k), around one of the ordering wave vectors K, or, 
alternatively, the use of Eq. (p^). We remind that the cor- 
relation length of the infinite system (above ^bkt) is well 
defined by Eq. (pi]); for a finite system, the same is still 
true for temperatures and lattice sizes large enough that 
the finite-size effects are negligible; otherwise Eq. (jl^) 
can be considered an ad hoc definition of ^. Provided 
that ^ < L/Q, the values of correlation length we got by 
the two methods were the same within the uncertainties. 
When finite-size effects become relevant, the results dif- 
fer from each other of about 4-8%; the first method being 
less reliable since Eq. ( p^ is valid for ^~^k^ ^ 1, a con- 
dition which cannot be satisfied beyond the wave vectors 
closest to K, as it can be seen in Fig. ||: the values of ^ 
reported in Fig. ^ for A = 0, are those obtained in the 
second way. The correlation length below ieKT, as well 
as in a neighbourhood above it, has been evaluated just 
to check where finite-size effects become relevant and in 
order to verify the finite-size scaling law ^ (x L. 

Fig. |o| shows the in-plane (x) and out-of-plane (x^^) 
sublattice susceptibilities, defined in Eq. (|^), as functions 
of temperature, for different values of A and of the lattice 
sizes. As expected, the susceptibility displays finite size- 
effects like the correlation length and the rule of thumb 
^ ^ L/6, for neglecting such effects also applies. The be- 
haviour of the out-of-plane sublattice susceptibility x^^j 
shown in the same figure, is also an interesting feature of 
the BKT transilipn, also present in the XXZ model on 
bipartite latticetS, which, of course, has no counterpart 
in the XY models. As expected, the absolute magnitude 
of x^^ increases, as lambda increases, the system becom- 
ing more isotropic and the out-of-plane fluctuations be- 
coming easier. However, for every value of A ^ 1, the 
easy-plane character of the system prevails at low tem- 
peratures and for t ^ Q, x^^ ^ 0; on the other hand, 
in the opposite limit, the effects of the anisotropy of the 
interactions disappear, all spins can fluctuate indepen- 
dently from each other and both the in-plane and out- 
of-plane susceptibilities approach the common value 1/3. 
Starting from the high temperature limit, and coming 
to lower temperatures, x increases whatever the value of 



A 7^ 1 is, ending to diverge at tsKT; of course, the smaller 
is the anisotropy, the longer x^^ follows the behaviour of 
X- The result is that, as A ^ 1, x^^ develops a sharper 
and sharper maximum. 
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B. Critical temperatures 

1. Order- disorder transition 

As we have seen in the previous sections the behaviour 
of the specific heat, chirahty and its susceptibihty as 
functions of temperature indicates the presence of an 
Ising-hke phase transition connected with the loss of the 
chirahty order for every value of A < 1 considered. Those 
thermodynamic quantities allow us to estimate imme- 
diately, at least approximately, the corresponding criti- 
cal temperatures, tc, since the transition appears rather 
sharp for all the values of A. A finite-size scaling analy- 
sis of these data shows that the features of such phase- 
transition are consistent with two-dimensional Ising ex- 
ponents. 
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FIG. 11. Finite-size scaling analysis for the staggered chi- 
rality k, for A = 0. Ising exponents /3 = 1/8 and v — 1 are 
used. Different symbols refer to the simulation box sizes as 
in Fig. |. 

In Fig. |ll| kL/^/" is reported, for A = and = 0.412, 
as a function of the reduced variable x —j:L}/'^; where r 
is l — t/tc for t < tc and 1 — tc/t otherwiseOO According 
to the scaling-hypothesis, close to the infinite lattice tc 
the order parameter k is given by 



(18) 



where, since for L oo the power-law singularities are 
to be reproduced, the limiting form of the function f{x) 
for a; — > oo and t < tc is 



f{x) cx X' 







(19) 



Therefore, reporting as in Fig. ID kL'^/" versus x. the 
data for the various lattice sizes do collapse onto a single 
curve, which is of course f{x) if the values of ^cj /3, and 
V are correct. In the present case using the critical ex- 
ponents of the two-dimensional Ising model, v = 1 and 
/3 = 1/8, a good agreement with the finite-size scaling 
law is obtained. Below the large- a; behaviour repro- 
duces the correct critical behaviour, which, in a double- 
logarithmic scale, is represented by a line with slope 



/3 — 1/8. Above tc the asymptoticj-bjehaviour shows the 
1/L decay of the order parameterc3 as i — > oo which 
means, according to Eq. (|l^), f{x) cx a;'^"^. This cor- 
responds to the line with slope —7/8 in the figure. In 
this way it has been possible to give the estimates of the 
critical temperature reported in Table |. 



2. BKT transition 

In order to estimate the critical temperature associ- 
ated with the BKT transition we relied on two methods: 
firstly, the fit of the correlation length and the in-plane 
susceptibility with Eqs. ( |6| ) and (p7|), using the MC data 
that are representative for their thermodynamic limit, 
i.e., for t ^ ^bkt; secondly, we used the finite-size scaling 
law 



X cx L 



2-r,(t) 



(20) 



which is valid for t ^ ^ekt, to obtain the scaling expo- 
nent r]{t)^ that satisfies rjit-^KT) — 1/4. Noteworthy, the 
latter method makes use of a different and independent 
data set. 

As for the first method, we have tested for BKT be- 
haviour, i.e., with Eqs. (|6|) and (p^, both the correlation 
length ^(i) and the susceptibility x(0- said above, in 
doing this only the reliable estimates of the thermody- 
namic limit of these quantities were kept, after the crite- 
rion ^{t) < 6L, in order to discard the data affected by 
finite-size effects. Of course, this has prevented us from 
obtaining useful estimates ^(t) and %(*) for t close to the 
transition temperature, due to the large computer-time 
required to simulate large systems. Indeed, when sim- 
ulating a larger lattice, besides the increase of the time 
needed for any move {oc N = L'^), we must also face the 
increase of MC fiuctuations and correlation time, so that 
much more sample configurations should be generated in 
order to keep uncertainties at a reasonable level. This 
would imply resorting to large-scale simulations which is 
well beyond our purposes. The results for the BKT tran- 
sition temperatures from the fits of the representative 
data are summarized in Table |l|; since x(<) is a direct 
outcome of the MC simulation, its values and the con- 
sequently fitted values of ^ekt are more accurate than 
those for ^{t), which must ideed be derived by fitting the 
MC outcomes for the k-dependent susceptibility. The 
uncertainties account both for the statistical error and 
for the instability of the fit against exclusion of the data 
points at the lowest temperature, where ^(t) L/6. 

As for the use of Eq. (20), in actual numerical calcu- 
lations it holds also slightly above ^ekt, when L is still 
smaller than the thermodynamic ^ and the system is al- 
ready correlated. In fact this scaling relation allows us to 
give an estimate of the parameter r]{t) and, by looking at 
which temperature such quantity attains the value 1/4, 
to have an independent check of the estimated critical 
temperature. In Fig. |l^, x/L'^^^ is plotted on a doubly 
logarithmic scale as a function of the lattice size L for 
A — 0.50. The data fall clearly on a straight line for 
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t < 0.40, the slope ol the lines being the corresponding 
values of 1/4 — 77; they already depart from linearity, in- 
stead, for t = 0.41. The values for the quantity 1/4 — 77 
we obtain by fitting the susceptibility data with Eq. ( |20| ) 
are reported, for the various values of A considered, in 
Table [J] . By interpolation, the values of ^bkt appear- 
ing in the fourth column can be computed. These data 
agree reasonably well with those obtained by fitting ^ 
and X with Eqs. and Eqs. (|7|), shown in Table ||; 
the trend to a slight overstimation of Ibkt for the higher 
values of A was already observed for the square-lattice 
case in Ref. nsl 
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FIG. 12. In-plane susceptibility over L'^''* vs L, for A = 0.5, 
at different temperatures around the critical one: t = 0.30 
(squares), 0.38 (up triangles), 0.39 (down triangles), 0.40 (di- 
amonds), and 0.41 (circles). 

For every value of A the critical temperature ic is 
significantly higher than the BKT transition tempera- 
ture iBKT, although their difference is not much larger 
than the uncertainties. Also in this respect the situa- 
tion is to that of the XY TAF, where there is no agree- 
ment in the literature on whether a single or a double 
phase-transition takes place. Nevertheless, the fact that 
tc ^ ^BKT for all values of A in a systematic way, makes it 
unlikely that the difference could just be due to statistical 
errors. 



it is consistent with a BKT transition with respect to 
the in-plane correlation length and susceptibility. The 
value of both the critical temperatures decreases with 
anisotropy strength, as shown in Fig. |l3|; this is consis- 
tent with the fact that the critical behaviour observed 
is connected with the planar character of the system, 
and that both the chirality and the orientational quasi- 
ordering are disturbed by the out-of-plane fluctuations of 
the spins which, at a fixed temperature, increase with the 
value of A. For the same reason the phase transition in 
the XXO model takes place at a temperature (t ~ 0.403) 
which is sensibly lower than that observed in the XY 
model (t ~ 0.505), where the spins are confined in the 
xy plane. As for the question of whether a single or two 
phase transitions are occurring, our results for the transi- 
tion temperatures (Fig. |3|) support the second hypotesis, 
consistently also with the most recent high-p«ecision MC 
simulations of the fully frustrated XY modelQ, where the 
existence of a new universality class is ruled out. 




FIG. 13. Transition temperatures for the Ising-Iike transi- 
tion (full diamonds) and for the BKT one (squares), the latter 
estimated through the BKT fit of the susceptibility, vs. the 
anisotropy. 



V. CONCLUDING REMARKS 



We have performed Monte Carlo simulations of the 
two-dimensional XXZ model on a triangular lattice at 
different values of the easy-plane anisotropy constant A. 
For every value of A considered, the situation appears 
quite similar to that observed in the XY triangular anti- 
ferromagnet, where frustration induces an order-disorder 
transition, associated with the two-fold additional degen- 
eracy of the ground state, and a BKT transition con- 
nected with the sublattice in-plane orientational order- 
ing. The critical behaviour turns out to be consistent 
with an Ising transition, for the internal energy, specific 
heat, chirality and the associated susceptibility, while 
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TABLE L Critical temperature tc associated to the chiral- 
ity phase-transition. 





A 






0.00 
0.50 
0.90 
0.99 


412 + 005 
0.400 ± 0.005 
0.355 ± 0.005 
0.320 ± 0.005 


TABLE IL BKT transition temperature as obtained fitting 
the in-plane correlation length £(t) and the in-plane suscepti- 
bility X against Eqs. (^) and (|l^). 


A 


tBKT (C fit) 


^BKT (X fit) 


0.00 
0.50 
0.90 
0.99 


0.402 ± 0.002 
0.391 ± 0.002 
0.345 ± 0.006 
0.306 ± 0.008 


0.403 ± 0.001 
0.388 ± 0.003 
0.344 ± 0.002 
0.305 ± 0.009 



TABLE in. Scaling exponent r^it) as obtained by direct 
finite-size scaling analysis of the in-plane susceptibility data 
according to Eq. (|20|). The values of tsKT given in the fourth 
column are obtained by interpolation of the function rjit) at 
the point r; = 1/4. 



A 


t 


1/4 -r, 




0.00 


0.300 


0.165 ± 


0.002 




0.400 


0.054 ± 


0.005 




0.405 


0.051 ± 


0.007 




0.410 


-0.06 ± 


0.01 


0.50 


0.300 


0.1654 ± 


0.0007 




0.380 


0.09 ± 


0.04 




0.390 


0.052 ± 


0.008 




0.400 


-0.25 ± 


0.03 


0.90 


0.300 


0.1486 ± 


0.0005 




0.320 


0.13 ± 


0.06 




0.340 


0.068 ± 


0.008 




0.350 


0.00 ± 


0.01 


0.99 


0.200 


0.204 ± 


0.009 




0.300 


0.131 ± 


0.02 




0.310 


0.08 ± 


0.01 




0.315 


-0.02 ± 


0.01 



0.407 ± 0.003 



0.391 ± 0.005 



0.350 ± 0.005 



0.314 ± 0.002 
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